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ABSTRACT 

^\ . The transient emission produced behind internal shocks that are driven by overtaking 

0^ I collisions of a magnetized, relativistic outflow is considered. A self-consistent model capable 

of describing the structure and dynamics of the shocks and the time evolution of the pair-and 
^ , gamma-ray distribution functions is developed and applied to gamma-ray flares in blazars, 

in the case in which gamma-ray production is dominated by inverse Compton scattering of 
external radiation (ERC). The dependence of the flare properties on magnetic field dissipation 
rate, intensity of ambient radiation, and the thickness of expelled fluid slabs is analyzed. It is 
shown that i) the type of gamma-ray flare produced by the model is determined by the ratio of 
the thickness of ejected fluid slab and the gradient length scale of ambient radiation intensity, 
\^ , ii) the radiative efficiency depends sensitively on the opacity contributed by the background 

' radiation, owing to a radiative feedback, and is typically very high for parameters characteristic 

to the powerful blazars, and iii) the emitted flux is strongly suppressed at energies for which the 
00 I pair-production optical depth is initially larger than unity; the time lag and flare duration in this 

energy range increase with increasing gamma-ray energy. At lower energies, flaring at different 
gamma-ray bands occurs roughly simultaneously, but with possibly different amplitudes. Some 
observational consequences are discussed. 
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Introduction 



The discovery of strong, variable gamma-ray sources identified with blazars by CGRO has motivated 
' many theoretical investigations concerning the emission from relativistic jets, and several models of 

^ i gamma-ray blazars have recently been developed (e.g.. Burns & Lovelace 1982; Dermer & Schlickeiser 

1993; Bloom & Marscher 1993; Manneheim 1993; Sikora, et al. 1994; Blandford & Levinson 1995 [BL95]; 
Ghisellini & Madau 1996). However, most of these efforts have been devoted to examine spectral properties 
of blazars using steady-state models, which are most suitable for exploring quiescent states. Although some 
progress in identifying the emission mechanisms in individual sources (e.g., Sambruna et al. 1997; Sikora 
et al. 1997), and constraining the structure and dynamics of relativistic jets on small scales (Levinson 
1996) has been made, it has became widely recognized that quantitative analysis of the growing body of 
variability data is crucial for advancing our understanding of the nature of relativistic jets further. Recent 
observational efforts (e.g., Reich et al. 1993; Maraschi et al. 1994; Wagner 1996; Buckley et al. 1996; 
Takahashi et al. 1996 ; Mattox et al. 1997; Wehrle et al. 1997; for a recent review see Ulrich 1997) to 
characterize the transient emission in blazars motivate detailed analysis of various variability mechanisms, 
and the development of time dependent models. 

Different processes may lead to time variability of blazar emission, including sudden changes in particle 
injection rate and/or magnetic field, changes in the bulk speed, and temporal variations of the intensity 
of background radiation in ERC models. This paper investigates the possibility that gamma-ray flares 
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observed in blazars are produced by internal shocks propagating in a magnetized, relativistic jet. The work 
presented in this paper extends an earher work by Romanova & Lovelace (1997; hereafter RL97). An 
outline of dissipative fronts is given in §2. In §3 we present the model and derive the basic equations. The 
results are described in §4. We conclude in §5. 

2. Dissipative fronts in magnetized, relativistic outflow^s 

Temporal fluctuations in the parameters of a MHD outflow lead to excitation of waves that steepen 
into shocks at some distance from the fluid ejection point, and the ultimate formation of a front consisting 
of a pair of shocks and a contact discontinuity across which the total pressure (kinetic plus magnetic) is 
conserved. The front propagates at a speed intermediate between that of the colliding streams and expands 
at a rate proportional to the relative velocity between the two shocks. Consequently, there is a net energy 
flow into the front which is balanced by either adiabatic cooling, owing to the front expansion, or radiative 
cooling, depending upon the conditions in the front. In the case of Poynting flux dominated outflows, rapid 
magnetic filed dissipation is required in order for a significant fraction of the outflow energy to be converted 
into kinetic energy behind the shocks. Efficient magnetic field dissipation may be achieved in strongly 
turbulent MHD flows (Thompson 1997) (an ordered magnetic field component and an associated Poynting 
flux can still be defined), or in the case of a sudden reversal of magnetic field lines (Romanova & Lovelace 
1992). 

We envision that the Lorentz factor of a fluid expelled from the central engine increases suddenly from 
r+ to r_ > r+, and denote by Tacc the characteristic time change of the outflow parameters (of order the 
dynamical time in the injection region). The ejection of the fast fluid is assumed to persist for time Tinj, 
after which it abruptly terminates, and the outflow Lorentz factor drops to lower values. As described 
above, this would lead to the creation of forward and reverse shocks that will propagate across the slow and 
fast fluid slabs, respectively. (The drop in Lorentz factor after time Tinj will lead also to the formation of a 
rare- faction wave behind the rapid slab.) The position To at which the front is created has been calculated 
by Levinson & van Putten (1997) using simple-wave analysis. As they have shown, in the limit in which the 
disturbance speed with respect to the rest frame of the boundary at which the fluid is injected (henceforth 
referred to as injection frame; this might be e.g., the rest frame of the central engine ) is highly relativistic, 
To = KCTacc^\T\, where Ta is the Lorentz factor associated with the Alfven speed with respect to the fluid 
rest frame, and k is a numerical factor that depends on the details of fluid injection (k = 1/3 when the 
4- velocity of the expelled fluid increases according to u = sinh[A + t/r]; t > 0). In the case of extragalactic 
jets, r+ is typically of order a few. Ta may range from ~1, in the case of a weakly magnetized flow, to 
> 10 beneath the annihilation radius in Poynting flux jets (BL95; Levinson 1996). Thus, if we associate 
cTacc with the gravitational radius of the putative black hole, we anticipate rg to lie in the range between 
10^^ and 10^* cm in blazars. 

The creation of fronts in non-relativistic, hydrodynamic flows has been considered by Raga et al. (1990) 
in the context of HH objects. The basic idea has later been generalized and applied to Poynting fiux jets in 
blazars by RL97. These authors proposed a scenario in which the gamma-ray and synchrotron outbursts 
often seen in blazars are ascribed to cooling of relativistic electrons accelerated inside highly dissipative 
fronts propagating in a magnetically dominated jet. The physical processes included in their model are 
inverse Compton scattering of external radiation, synchrotron and SSC emission. In their treatment, 
however, the spectrum of emitting particles is assumed a priori, and is characterized by two break energies 
which are allowed to evolve with time. The structure of the front is not calculated self-consistently, but 
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rather the front is assumed to expand freely with the speed of sound. As we show, the front structure 
can be significantly altered by radiative losses. Moreover, their model does not account for absorption of 
gamma-rays escaping the front by pair production on ambient photons ahead of the front. As shown below, 
this process may considerably affect the emitted spectrum, particularly in the powerful blazars. 

In this paper we generalize RL97 treatment in several important ways. Firstly, the time evolution of 

the energy distribution of pairs and gamma-rays in the front is calculated sclf-consistcntly, by numerically 
integrating the coupled kinetic and MHD equations. Secondly, the c^ffcct of radiative losses on the evolution 
of the front structure is taken into account. Thirdly, pair cascades inside and ahead of the front are 
accounted for properly. Fourthly, finite length outbursts are considered. And finally, we use different initial 
conditions. Our model involves the following kc;y parameters: i) the dissipation rate of magnetic energy, ii) 
the maximum injection energy of electrons, denoted by Eemaxi hi) the fraction rj of dissipation energy that 
is tapped to the injection of electrons to Egmaxt and iv) the ratio of the thickness of ejected fluid slab and 
the gradient length scale of the intensity of background radiation. We suppose that on scales of interest the 
outflow material is dominated by e^ plasma, and ignore possible effects associated with entrainment of ions 
(which are likely to be present in the MHD wind confining the jet) on the front dynamics. In the present 
work, we consider only gamma-ray flares produced by the ERC mechanism. The extension of this work to 
encompass synchrotron and SSC flares will be presented in a follow-up paper. 



3. The model and basic equations 
3.1. Front dynamics 

The stress-energy tensor of a magnetized flow can be written as a sum of two contributions: the 
energy-momentum tensor associated with the electromagnetic field, T^, and the energy- momentum tensor 
associated with the system of particles, T^". T^" and the particle fiux. A'''*, are related to the electron 
distribution function, defined in §4.2 below through, 

N^ = nu^ = Jp'^f,{p,x)^, (1) 

and 

Tr = Jp'p-fsiP,x)^, (2) 

where n is the proper density and u" the fluid 4-velocity. In terms of the proper gas pressure, p, specific 
enthalpy, h, magnetic induction, 6" = F^^u^, with F"'^ being the dual electromagnetic tensor, and rest 
mass density p = men, the stress-energy tensor takes the form 

T^'' = T,^"" + TZ = ph*u''u''' + p*g'''' - b^'b" . (3) 

Here gap denotes the metric tensor, h* ~ h + fe^/ (pc^), p* = p + &^/2, and b^ — 6"6q. 

Let Q and S'^ denote the rate of change of the electron density due to e^ pair creation and annihilation, 
and the rate of change of the u component of the stress-energy tensor resulting from radiative losses, 
respectively, and denotes the source term associated with magnetic field dissipation. The fluid equations 
can then be expressed as, 
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5aAf" = Q, (4) 

The terms Q and are expressed in §4.3 below in terms of the Boltzmann coUision operators. The 
evaluation of these source terms at any given time involves the integration of the kinetic equations (eqs. 
10 and jl^) for the pairs and photons. 

We consider the evolution of a front produced by a collision of fast and slow fluids having velocities f3^ 
and and corresponding Lorentz factors r_ and r+, respectively. We suppose that each shock decays 
after it crosses the corresponding fluid slab, at which point energy deposition into that slab ceases. At times 
shorter than the shock crossing times, the set of equations (||) can be integrated over a portion of space-time 
containing the front, using the generalized Gauss's theorem, to obtain a set of ordinary differential equations 
governing the time evolution of the front quantities. The derivation of the front equations is given in the 
appendix. In what follows, subscript minus (plus) refers to quantities leftward (rightward) of the contact 
discontinuity, and subscript / denotes quantities inside the front. For simplicity, we shall restrict our 
analysis to the case in which the colliding fluids have the same proper density, pressure, and magnetic 
field, viz., — p+ = 6+ = b—. In this case a symmetric front will be created, i.e., nj+ — 
Tj'^ = Tj!''^. The problem is then characterized by seven independent variables: Tj°, T^^, Nf, Bj, f3s+, Ps-, 
Pc, where (3s± are the corresponding shock velocities, and Pc is the front velocity (i.e., the velocity of the 
contact discontinuity). Since the proper energy density of the fast and slow fluids is the same, the energy 
flux of the fast outflow, as measured in the injection frame, is larger than that of the slow outflow by roughly 
a factor of (r_/r_|_)^. Thus, this case corresponds to a situation wherein the central engine undergoes an 
outburst during which energy is being injected into the system. An alternative possibility, which will not be 
considered here, is that the fluid parameters change in such a way as to sustain the associated energy flux 
unchanged. This would lead to the creation of an asymmetric front, which in some circumstances can give 
rise to a somewhat higher radiative efficiency. To simplify the analysis further, we suppose that the proper 
density, pressure and magnetic field in the front are homogeneous, that d+ = d- — d, and that the magnetic 
field dissipation rate is proportional to the injection frame magnetic field, -B/, and inversely proportional to 
the front length AX, viz., Qi,A.X = —ajjBf, where ab is a constant. The latter assumption is not crucial, 
but is convenient since it admits steady-state solutions in the adiabatic case (i.e., in the absence of radiative 
losses) which simplifies the choice of initial conditions (cf §3.5). As we have verified, the essential features of 
the solution are insensitive to the choice of the form of the term associated with magnetic field dissipation. 



With the above simplifications eqs. (A7) reduce to 



AX^iNf) = AXQ - Nf{p,+ ~ps^)+ N+iP,+ - /3+) - A^_(/3,„ - /3_), (5) 

AX^iBf) = -a^Bf ~ BfiPs+ - Ps-) + B+{P,+ -/?+)- i?-(/3s- - P-). (6) 

AX^Tf - AX5° - r™(/3,+ -ps-) + {T°"Ps+ - Tf ) - (T^^Ps- - T^"), (7) 

AX^Tf = - Tf{p,+ ~ps-) + {Tl^Ps+ - T^l ~ (T^^A- - T^"), (8) 

where x° = ct, t being the time measured in the injection frame, and where the shock and front velocities 
obey the algebraic equations, 

(r;° - T°°)/3,+ + (r;° - r°°)/3,_ = 2t;^ - - r°^, (9) 
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(t;^ _ T°^)/3^^ + (t;^ _ T°^)/3,_ = 2t;^ - r_^^ - t^^, (lo) 
(r - i)t;"/32 - [r(r;° + b)/2) - (r - i)(r/ + b))](i,, + t;" = o. (ii) 

The last two terms on the RHS of equations (^) account for the flux of the corresponding quantity 
incident into the front through the shock surfaces, whereas the second term on the RHS describes the rate 
of change of the corresponding front quantity due to the front expansion. Finally, the expansion of the front 
is governed by the equation, 

= (12) 
and the crossing times of the forward and reverse shocks, denoted by tj^±, are given implicitly by 



R± 



d±/c = ± / {I3,± - fi±)dt. (13) 
Jo 

As discussed earlier, the above equations are valid only for times shorter than min {tR+,tR-} ■ The equations 
describing the evolution of the system at later times can be derived in a similar manner. Formally, when 
t > tjij^ we replace I3s+ by Pc in eqs. (|5|)-(|ll|) , and set A^+, _B+, and T^"^, to zero. Likewise, when t > t^- 
we set /3s- — (3c, = B- = T^^ = 0. In doing so we ignore the rare- faction waves that will be produced 
at the edges of the slabs due to the pressure gradient there. 



3.2. Kinetic equations 

Let fi,{pe,x) and f^{p^,x) denote the distribution functions of electrons (we do not distinguish 
electrons from positrons, designating both by subscript e) and gamma-rays as measured in the injection 
frame, respectively. The corresponding Boltzmann equations can be written as 

p':dJe^Cc{f,P,x)+anj, (14) 

and 

pi;9^^ = C^(/,p,a;). (15) 

The operators Ce and on the RHS of the above equations represent the interaction between pairs and 
gamma-rays, and are given explicitly in BL95. They are restricted to the condition 

(Ce + Cj)d^p = 0, (16) 

by virtue of energy conservation. The injection operator, Cinj , represents the change in energy and 
momentum of electrons (positrons) due to their interaction with the large scale electromagnetic field; that 
is, as a result of shock acceleration, magnetic reconnection, or stochastic acceleration owing to absorption 
of nonlinear plasma waves. Since the injection operator preserves the total electron number density, it must 
satisfy, 

a„,^ = o. (17) 



The angular distribution of pairs and gamma-rays is expected to be strongly beamed along the direction 
of propagation of the front in the injection frame. Thus, we can approximate the distribution functions 
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as: fa ^ Pa'^na{Ea, x)S{^a); a — (e,7), with /i^ being the cosine of the angle between the momentum of 
a species a and the flow velocity, and na{Ea) the corresponding number density of a species a per unit 
energy. Integrating eqs. (^) and ( |l5|) over a region encompassing the front in the manner described in the 
appendix, assuming again that the distribution functions are homogeneous inside the front, and averaging 
over /ict yields to a good approximation, 



d 

AX—Uef - I {Ce+ C^no)dx ~ ne/(A+ - A-) + rie+(/3,+ - 1) - n^- - 1), (18) 
dx 



d 

AX—n^f ^ I C^dx - n^f{l3s+ - 13s-) + n^+{l3s+ - 1) - n^-{l^s- ~ 1). (19) 



Now, the electron density in the fluids exterior to the front, viz., ne+ and ne— , are given as input, 
and can be absorbed into the definition of the injection operator. For convenience we redefine the injection 
operator as follows: / Cmjdx J C'mjdx — J Cinjdx + ne+(/?s+ — 1) — ne-{l3s- — 1). By employing eq. 
( p^ ) we arrive at, 

' dx I dE, = N+ (/3,+ -/?+)- N- ( A_ ~ /3_ ) . (20) 



This condition simply reflects the fact that particle acceleration preserves the total number density of 
particles incident into the front. A second constraint on Cmj is obtained by multiplying eq. ( p^ ) by i?e, 
and then integrating over dEe- 

dx J E.anjdE, = T°;(/?3+ + j - / (21) 

where T°j is given by eq. (||). Note that the RHS of the last equation equals the fraction of total energy 
flux dissipated inside the front. 

We must also determine and n^_. Under the approximation of perfect beaming invoked above 
= and n~f^ — n-yf, since the only source of beamed radiation is the front itself. With the above 



results eq. ( |19| ) simphfies to, 

-^n,f = + iAX)-^n,j{(3s- ~ 1). (22) 

Again, we assume that for t > tj^± the time change of the distribution functions is dictated by eqs. ( p^ 
and ( |2^ ) with and /3s+ replaced by /3c and n^± set to zero. 

The emergent gamma-ray spectrum depends on the pair production opacity contributed by the ambient 
radiation field ahead of the front. Gamma-rays for which the pair production optical depth to infinity 
largely exceeds unity, will be converted into lower energy pairs and gamma-ray via pair cascades upstream 
the forward shock. (We note that the escape probability of gamma-rays at a given energy from the front 
may be considerably different than the probability that they will escape the system to infinity. In fact, 
the pair production optical depth of the front itself will be smaller than unity as long as its axial length is 
smaller than the corresponding gamma-spheric radius.) This leads, as shown below, to a strong suppression 
of the emitted flux at the corresponding energies. The calculation of the emitted spectrum proceeds as 
follows: At every time step we integrate the equations describing the spatial evolution of the cascade, 

^n,(£;,) = C„ (23) 

-—n,{Ee) = Ce, (24) 
a Inr 
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starting at the current position of the forward shock, r(t) = J'o + / Ps+dt, and subject to the boundary 
conditions ne[r(i),£'e] = ne+^ n~f[r(t), E^] = n^f(t,E^). Here n~f{r,E^) is the number density per unit 
energy of escaping gamma-rays ahead of the forward shock. The assumption underhne these calculations 
is that the cascade develops on a sufficiently short time scale to render any retardation effects negligible. 
We also ignore possible alterations of the front structure due to momentum exchange between emitted 
gamma-rays and the slow outflow (i.e., upstream of the forward shock). The inclusion of this effect 
complicates the numerics substantially and is beyond the scope of our analysis. We anticipate, though, that 
such alterations of the front structure will not affect the emission characteristics significantly. 



3.3. Determination of the source functions 

The source term describing the rate of change of the electron density in the front due to formation of 
pair cascades, is determined by integrating eq. (|T^) over energy, and using eqs. (|o|) and One then 
obtains, 

Q = 2y CedEe ^2 J Kpp{E^)7i^fdE^, (25) 

where Kpp{x, Ej) is the pair production opacity, and the second equality has been obtained using eq. (4.12) 
of BL95. Likewise, the radiative energy loss rate is obtained by taking the first moment of eq. (|l8|), i.e., 
multiplying the equation by E^ and then integrating over dEg, and by using eqs. ( pT| ) and (^): 

S° = J EeCedEe. (26) 

Under the assumption of perfect beaming the momentum loss rate equals the energy loss rate, viz., 

— cS° . By employing eqs. ( |l6| ) and (|2|) we can rewrite the rates associated with radiative losses in the 
form 

QrpOO 

S-/c = S° = -—^ - (AX)-iT°°(l - /3,_), (27) 
ox° ' 

where T°° = J E^n^dlnE^ is the gamma-ray energy density. Note that the latter equation is essentially 
the energy equation for the system of gamma-rays. 



3.4. The injection operator 

As stated above, the injection operator represents the energy redistribution of electrons (positrons) 
resulting from their interaction with the electromagnetic field, e.g., due to Fermi acceleration, magnetic 
reconnection, or absorption of nonlinear plasma waves excited by the shock. The physics of these processes 
is yet to be understood better before Cinj can be derived from first principles. Here we settle for a simple 
prescription in which a small fraction of the pairs incident into the front are injected to some maximum 
energy, assumed to be fixed in the observer frame, and the rest are redistributed at much lower energies. 
The maximum injection energy, Eemax, and the fraction 77 of the dissipation energy that is being injected 
to Eemax are treated as free parameters. To be concrete, we choose 

C,nj = niS{E^ - Eo) + n25{E^ - E^^^x) (28) 

with rii, and Eq determined from eqs. ( pO| ) and (|2l|), and the requirement that n2Eemax comprises a 
fraction rj of the total dissipation energy. Other forms of Cinj can be readily derived, for example, a thermal 
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distribution with a power law tail. We find though that the results are insensitive to the exact form of the 
injection operator provided that particle acceleration is sufficiently efficient, in the sense that a significant 
fraction of the injection energy is carried by electrons having energies near E^max- K should be noted that 
the assumption that E^max is fixed in the observer frame is unrealistic. A more realistic prescription would 
be to take E^max to be fixed in the rest frame of the front. Unfortunately, this complicates the numerics 
considerably. However, variations of E^max due to the changing front velocity should not affect significantly 
the evolution of the gamma-ray distribution function at sufficiently low energies. 



3.5. Initial condition 

In order to integrate the above equations, one must specify the initial position and structure of the 
front, and the initial energy distribution of electrons. A simple and convenient choice would be to use 
the structure of an adiabatic front. Below, we employ the front structure calculated by Levinson & Van 
Putten (1997). Specifically, we first solve algebraically eqs. (|5|)-(|lT|) in the absence of radiative losses (i.e., 
with Q — S° — ^ ah ^ 0), and assuming steady state {d/dx° = 0), for the chosen input parameters. 
The front quantities thereby obtained then serve as initial values for the integration of the radiative front 
equations. This choice is viable if the creation of the front occurs on a timescale much shorter than the 
radiative cooling time. An initial electron distribution subject to the restrictions that the number density 
and total energy (temperature) equal those computed for the adiabatic front is also specified. The initial 
thickness of the front, AX(t = 0), is taken to be small enough, so that the energy loss rate due to cooling 
of the initial population of electrons is well below the rate of energy deposition in the front. This renders 
the results highly insensitive to the choice of initial electron distribution. 



4. Results 

Equations (|)-(©, (©, (||)- (||) and (H) have bee n integrated numerically using the initial 
conditions discussed in §3.5. The source functions Q and 5'° have been computed using eqs. ( ^5| ) and (|2^). 
We have verified that the solution is indeed highly insensitive to the choice initial electron distribution 
when the initial length of the front is sufficiently small, as stated in §3.5. The equations have been modified 
in a manner described in §3.1 after each shock crossing (first at i = min {ta+TtR-} and again at i = max 

t/f_}), and the integration continued using the modified equations. In the following examples, the 
Lorentz factors of the slow and fast fluids and the rest frame Alfven 4-velocity have been taken to be 5, 20, 
and 10, respectively, and the initial electron distribution has been taken to be Maxwellian with appropriate 
temperature and density. A rapid magnetic fleld dissipation has been invoked with af, — 0.5. The standard 
soft photon intensity (consisting of a broken power law with a steeper slope below about 0.5 KeV) defined in 
Levinson & Blandford (1995) has been adopted for the calculations. To be concrete, we used the following 
form for the intensity of external radiation: 

h{Es,t) ^ ^^g{Es); lO-"^ < < 0.1, 

with r{t) =ro + J Pcdx°, and g{Es) cx {Es/Eo)-^^^[l + {Es/Eo)-^]; J g{Es)dEs = 1. Here Eo is the break 
energy and is taken to be 0.5 keV, and e << 1 is the fraction of nuclear luminosity that is reprocessed or 
scattered by surrounding gas. As a check, we integrated eq. ( p7| ) to obtain T°° and compared the result 
with total flux emitted (i.e., J E^n^dlnE^) at every time step. The agreement was typically better than 
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2%. 

To illustrate the dynamics of the system we first consider the limit of continuous ejection; that is, 
d/ro — > oo. The time evolution of the front quantities is depicted in fig. 1, where the 4-velocity of the 
contact discontinuity and the two shocks (upper left panel), the proper density (upper right panel), the 
total pressure (bottom left panel), and the proper magnetic pressure (bottom right panel) are plotted 
against \og{ct/ro), t being the injection frame time, for d/vo — 100, (eLs)45/''i6 = 1 (solid lines) and 
(eLs)45/'"i6 = 10~^ (dashed lines). Here (6^^)45 is the scattered luminosity in units of 10''^ erg s~^, and 
Tie, = ro/(10^^cm). The corresponding radiative efficiency, defined as the fraction of dissipation power 
radiated by the front, and given explicitly by, 

^^^^ ^ Jn.,iE^:t){l- P,+)dE^ 

J dx J EeCinjdEe 

is exhibited in fig. 2. As seen from figs. 1 and 2, the front velocity and expansion rate, initially equal those 
of an adiabatic front, decrease as the emitted energy flux rises. The rest mass density and total pressure 
increase correspondingly. After the peak emission is reached, the front starts accelerating and adiabatic 
cooling becomes gradually more important until, ultimately, the initial structure and velocity of the front 
are restored. The increase in number density is partly due to pair production and partly due to enhanced 
compression resulting from the drop in expansion rate. The change in total pressure is primarily due to 
magnetic field compression. Quite generally, we find that when t]i± > l/c, the maximum of the flux is 
reached at a distance I from the creation radius r^, where / is the gradient length scale of the intensity 
of ambient radiation at {I ^ in this example), and the duration of the flare corresponds to several 
times / (see below). When tR± < l/c the timescale of the flare is determined by the shock crossing times. 
The evolution of the system in this case is essentially the same as described above for times shorter than 
tR±- At later times, however, energy supply to the hot outflow terminates, and the outflow begins to cool 
radiatively and decelerate until radiative losses become small. The radiative efficiency, namely the ratio of 
radiated and dissipated energy fluxes increases with increasing values of (eis)45/''i6i as seen from fig. 2. It 
approaches 80 percent near the peak in this example for {eLs)i^/riQ = 1, and is generally around this value 
for typical parameters of gamma-ray blazars. 

The rate of energy dissipation in the front is presented in fig. 3, and it is seen that it increases as the 
radiated flux increases. The reason is that radiative losses lead to deceleration of the front and shocks (see 
fig. 1) and, consequently, enhancement in the rate of energy deposition behind the shocks, particularly 
the reverse shock, as can be inferred from eq. (|^). This positive feedback has important implications for 
the emitted spectrum, particularly in the presence of large pair production opacity, that are discussed 
below. Fig. 4 depicts typical light curves computed for d/ro > 1 and different values of (eLs)45/?'i6- In this 
figure the total apparent gamma-ray luminosity is plotted against the time measured by a distant observer, 
tabs = /o [1 ~ Ps+{t') cos 0\dt' , where [3s+ is the velocity of the forward shock, and is the angle to the line 
of sight (6* = in this example). The strong dependence of the flare intensity on the parameter (6X5)45 /rig 
is evident. This is one consequence of the positive feedback mentioned above. The decay of the flare in this 
regime depends on the radial variation of the intensity of background radiation and the expansion rate of 
the front, and is typically longer and more gradual than the rise. Further, the rise time increases slightly 
with increasing L^c^/riQ. The reason is that larger opacity results in smaller shock and front velocities 
during peak emission (see fig. 1) and, therefore, smaller beaming factor. From fig. 4. it is seen that the rise 
time is of order {I /2c)T~^{tpeak) — 10~^ro/c, where Ts+itpeak) is the Lorentz factor of the forward shock 
near maximum flux, and that the flare duration is several times longer. For Tq between 10^^ and 10^® cm 
(cf. §2) this corresponds to a rise time in the range between several minutes and several weeks. However, 
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gamma-rays of a given energy cannot escape from a radius smaller than their gamma-spheric radius, so 
that the duration of a flare in a given band is also limited by the pair production opacity. In general, we 
find that for parameters typical to the powerful blazars the flare duration in the EGRET band can range 
from several hours to several weeks for small viewing angles (i.e., 6 < T~^), when d/vo > 1. The shape of 
the flare is substantially altered when d/ro becomes sufficiently small. This is demonstrated in Fig. 5 where 
light curves computed for L45/ri6 = 1 and different values of d/ro are displayed. As seen, for values of 
d/ro smaller than unity the time scale of the flare, as measured by a distant observer, is of order d/c\ the 
light curves exhibit a steep decline at times tobs > d/c. For sufficiently thin slabs {d/ro < 2 x 10~^ in this 
example) the decay time of the flare is comparable to the rise time or even slightly shorter. 

In order to examine the relationship between the emergent fluxes in different high-energy bands 
predicted by our model, we divided the energy interval into several subintervals and followed the time 
evolution of the flux in each band. Quite generally we find that i) only a relatively small fraction of the 
radiated energy is emitted above the gamma-spheric energy at which the pair-production optical depth at 
r ^ To + l (the front position at peak emission) equals unity, ii) the time of peak emission and flare duration 
increase with increasing gamma-ray energy in this energy range, and iii) below that energy variations in 
the fluxes at different X-and gamma-ray bands occurs roughly simultaneously, but not necessarily with the 
same amplitude. This behavior is also a consequence of the positive feedback discussed above. An example 
is given in fig. 6, where the time evolution of the energy fluxes emitted in four equally spaced subintervals 
in the energy interval 5 MeV to 50 GeV, for d/l > 1 is presented. Similar behavior has been found in the 
case d/l « 1. This result is in contrast with the simultaneous flaring predicted by the model of RL97, 
which does not take into account gamma-ray attenuation by pair production on external photons. 

We also checked the dependence of the front dynamics and emission on the rate of magnetic fleld 
dissipation. We flnd that as Ua increases above unity, the front expands more rapidly, adiabatic cooling 
becomes more important and, hence, the radiative efficiency drops. The rapid expansion results in a slower 
decay of the flare, but otherwise the shape of the light curve is not altered significantly. 

5. Conclusions 

This paper considers the production of gamma-ray flares through inverse Compton scattering of 
external radiation by pairs accelerated behind internal shocks, that are driven by temporal fluctuations 
in the parameters of a magnetized, relativistic outflow. A self-consistent model capable of describing 
the dynamics of the front and the time evolution of the angle averaged pair and gamma-ray distribution 
functions, has been developed and employed to calculate gamma-ray light curves in blazars. 

The main results and conclusions of this study are: 

1) The shape and timescale of the flare depend on the ratio of the thickness of ejected fluid slab, d, and 
the gradient length scale of the background radiation intensity at the radius of shock formation, I. When 
d/l is sufiiciently small, such that the propagation time of the shocks across the fluid slab is shorter than 
l/c, the flare duration, as measured by an observer at small viewing angle, is ~ d/c, and the shape of the 
light curve is roughly symmetric. For larger values of d/l the flare reaches its maximum at a distance I 
from the creation radius. It then decays gradually until time d/c, after which the flux declines steeply. For 
d/l > 1 the shape of the flare is determined essentially by the radial variations of the intensity of ambient 
radiation and the expansion rate of the front; typically, the rise is fast compared with the decay. The above 
results suggest that different types of flares will be produced in different sources, or even in the same source 
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under different conditions. This expectation appears to be consistent with the variety of types of EGRET 
flares (e.g., Hartman et al. 1993; Kniffen et al. 1993; Mattox et al. 1997) observed in gamma-ray blazars. 
In the event of impulsive outflow ejection (i.e., ejection time of order the dynamical time, — Tacd 
cf. §2) the axial length of the fluid slab can be as small as the size of the central engine. If the latter is 
associated with the gravitational radius of the putative black hole, r^, then flare durations as short as 
Vg/c ~ 50[Afs_f//(10''MQ)] sec are plausible. Successive ejection of jets with short duty cycle may also lead 
to a substructure in the light curve of a longer duration outburst (as often seen in GRBs). In general we 
find that for parameters typical to extragalactic jets, namely F of order 10 and comoving Alfven 4-velocity 
between and 10, the flare timescale can lie in the range between several minutes and several weeks (for 
small viewing angles) . We note that the shape of the light curve may be significantly altered when the SSC 
process becomes important (cf. RL97). The study of SSC fiares is left for future investigation. 

2) The peak flux and flare intensity depend sensitively on the Thomson opacity near Tq, specifically, on 
the value of Lg/vo, where Lg is the fraction of nuclear luminosity that is reprocessed or scattered across the 
jet. This behavior is a consequence of a positive feedback that gives rise to a strong enhancement of the 
energy deposition rate during peak emission. For parameters typical to the powerful gamma-ray blazars, 
the radiative efficiency, namely the fraction of dissipation energy that is radiated by the front is typically 



3) The amplitude of variations of the flux emitted at energies for which the pair-production optical 
depth is initially larger than unity is much smaller than that at lower energies. The flare, in this energy 
interval, propagates from low to high energies and its duration increases with energy. At lower energies, 
flaring in different X-and gamma-ray bands occurs roughly simultaneously, but with possibly different 
amplitudes, depending on the spectra of ambient radiation and injected electrons. The claimed invariance 
of the gamma-ray spectrum during the strong EGRET flare observed in PKS1622-297 (Mattox et al. 1997) 
is not in conflict with that model prediction, since the sensitivity was insufficient to resolve any time lags, if 
present, between the fluxes in the two energy subintervals (below and above 300 MeV) analyzed by those 
authors, particularly in view of the small amplitudes anticipated at high EGRET energies. Despite the 
small amplitudes, it may be possible to observe time lags between the emission at hard (but still well below 
TeV where absorption by the IR background is strong) and soft gamma-ray energies with a next generation 
gamma-ray telescope. The detection of such lags is an important test for this model. One implication 
of the above results is that in some cases strong flares may be observed in some energy bands while the 
simultaneous variations in other bands may appear modest or even small. The simultaneous X-ray/TeV 
flare reported recently for Mrk 421 and the lack of signiflcant variations of the EGRET flux (Macomb, et 
al. 1995; Takahashi et al. 1996) is an example. The hard X-ray flux may be due to either inverse Compton 
emission of the thermal electrons or synchrotron emission of the highest energy electrons. Alternatively, this 
event may be explained as an SSC flare resulting from changes in the tail of the distribution of non-thermal 
electrons (Takahashi et al. 1996). 

I thank the referee, R.V.E. Lovelace, for useful comments. Support by Alon fellowship and a TAU 
Research Authority grant is acknowledged. 
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A. Appendix: Derivation of the front equations 



Equations (||) 



are of the form. 



d^A^ = B. 



(Al) 
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Integrating the above equation over a region in space-time enclosed by the world lines associated with the 
left shock, the contact discontinuity and the lines defined by the equations t — to\ t = to + At, where to is 
some fiducial time [see fig. (7)] , and using the generafized Gauss's theorem, 

M'dPY.^ = j Bd^n, (A2) 

where d'^il is a volume element and is the surface area of the 3d surface enclosing the region of 

integration yields, 

{A^/3,_ -A1+ A}_ {xc) ~ Aj_ {xc)f3c}cdt + / A}_ {to)dx 

o Jxa{to) 

- A°f_{to + At)dx = cdt B{t)dx, (A3) 

where it has been assumed that A^^ and B depend solely on the axial coordinate (i.e., along the 
direction of motion). In the limit At — ^ 0, and using the expansions Xc{to + dt) ~ Xc{to) + Pc{to)dt] 
x-{to + dt) ~ x^ito) + f3s-ito)dt, we arrive at, 

r^- dA°f z-^- 

A'L(3s~~Al-A}_{x-)f3s-+A}_{xc)+ B^dx = 0. (A4) 



Likewise, integrating eq. (Al) over a portion of space time containing the region rightward to the 
contact discontinuity yields, 

dA° 

-A;/3,+ + A^ + A^+(.t+)/3,+ -Af+(x,)+y —I±dx-J B+dx = 0. (A5) 

Under the assumption that Af± are homogeneous inside the front, (i.e., independent of the axial coordinate 
X, but may still have different values on each side of the contact discontinuity), we obtain the following 
differential equations for Aj:±{t) 

± {AX±)-\Alf3s± -Al- A'}^(3s± + A^±) - ^ + B± = 0, (A6) 

where AX±{t) — ± j^*"^^ dx is the distance between the forward (reverse) shock and the contact 
discontinuity. 



Substituting the front quantities for one obtains. 



(dNf± 



± [N}±{f3s± - Pc) - N±i(3s± - /3±)] + AX± - Q±j = 0, 

± [Bf±Ws± - Pc) - B±{f]s± - P±)] + AX± (^^^ - Qb±) = 0, (A7) 

dx° 



±[(tOO/3,± - rO|) - (Tf /3,± - Tf )] +AX±i—i±-Sl]=0, 



CgrpOx \ 

Here Bf± — bf±^c is the magnetic field as measured in the injection frame, Ps± are the corresponding shock 
velocities, and /3c is the velocity of the contact discontinuity. 
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Fig. 1. — Time evolution of front quantities. Shown are (a) the 4- velocity of the contact discontinuity 
(labeled Uc), forward shock (labeled Us+) and reverse shock (labeled U-), (b) the proper density, (c) the 
total pressure, and (d) magnetic pressure as a function of ln(ci/ro), for {eLs)45/rie = 1 (solid line) and 10~^ 
(dotted line). Here (ei/s)45 is the fraction of ambient luminosity scattered across the jet in units of 10^^ ergs 
s~^, and ri6 is the radius of shock formation in units of 10^^ cm. 




Fig. 2. — Radiative efficiency, defined as the fraction of dissipation power radiated by tlie front (see text), 
versus log of injection frame time. 
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Fig. 3. — Time evolution of the dissipation power. The dissipation flux is normalized to the bulk energy flux 
of the fast fluid. Curves are labeled by values of log[(eI/s)45/''i6] 



-17- 



LIGHT CURVES 
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Fig. 4. — Typical gamma-ray light curves produced by the model for d/rg > 1. Shown is the total apparent 
luminosity in units of Lj, the power of the ejected (fast) outflow, against ctobs/i~o, tots being the time 
measured by distant observer (see text for details). The numbers that label the curves are the values of 
(e^s)45/r-i6- 
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Fig. 5. — The same as fig. 4, but for (eLs)45/ri6 = 1 and different values of d/vo (label the curves). 
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Fig. 6. — Spectral evolution: displayed are the total apparent luminosities in four equally spaced energy 
intervals in the range 5 McV to 50 GeV (indicated in window a; energies arc given in units of GeV), 
for d/vo > 1 and (a) (el/s)45/n6 = 1, (b) (eLs)45/n6 = 10~\ (c) , (eLs)45/n6 = 10"^ and (d) 
(eL,)45/ri6 = 10-3. 
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Fig. 7. — Schematic diagram showing the world-Unes of the forward shock (/3+), reverse shock (/?-), and 
the contact discontinuity (/3c), and the corresponding contours of integration (dashed hnes). The region 
marked by the minus (plus) sign corresponds to the portion of the front leftward (rightward) of the contact 
discontinuity surface. 



